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We present a macroscopic calculation of coherent electro-magnetic radiation from air showers 
initiated by ultra-high energy cosmic rays, based on currents obtained from three-dimensional Monte 
Carlo simulations of air showers in a realistic geo-magnetic field. We discuss the importance of a 
correct treatment of the index of refraction in air, given by the law of Gladstone and Dale, which 
affects the pulses enormously for certain configurations, compared to a simplified treatment using a 
constant index. We predict in particular a geomagnetic Cherenkov radiation, which provides strong 
signals at high frequencies (GHz), for certain geometries together with "normal radiation" from the 
shower maximum, leading to a double peak structure in the frequency spectrum. We also provide 
some information about the numerical procedures referred to as EVA 1.0. 



I. INTRODUCTION 

The aim of our work is to provide a realistic calculation 
of radio emission from air showers, which might be used 
finally to analyze and understand the results from radio 
detection experiments (LOPES 0, Q, CODALEMA @], 
LOFAR Q), and the new set-ups at the Pierre Auger 
Observatory (MAXIMA @, AERA @|). 

There are two ways to compute the electric fields cre- 
ated by the moving charges of air showers: the "macro- 
scopic approach" adds up the elementary charges and 
currents to obtain a macroscopic description of the to- 
tal electric current in space and time, which is the source 
of the electric field obtained from solving Mawell's equa- 
tions. The "microscopic approach" computes the fields 
for each elementary charge, and adds then all the fields 
(with a large amount of cancellations). 

Already in the earliest works of 0413 > a macroscopic 
treatment of the radio emission was proposed, but at 
the time the assumptions about the currents were rather 
crude. There is recent progress concerning the macro- 
scopic approach. In 2007, we performed calculations al- 
lowing under simplifying conditions to obtain a simple 
analytic expression for the pulse shape, showing a clear 
relation between the pulse shape and the shower profile 
fill ]. This allows, for example, to determine from the ra- 
dio signal the chemical composition [l2| of the cosmic ray. 
The picture used was very similar to the one in Ref. [8|, 
which has been refined by using a more realistic shower 
profile and where we calculate the time-dependence of 
the pulse. Recently it was confirmed that the pulse pre- 
dicted in the microscopic description [13] agrees with 
the predictions following from the macroscopic picture as 
shown in [l5j . 

In Ref. [Ill , we advance further by computing first the 
four-current from a realistic Monte Carlo simulation (in 



the presence of a geo-magnetic field), and then solve the 
Maxwell equations to obtain the electric field, while con- 
sidering a realistic (variable) index of refraction, given by 
the Gladstone-Dale law as 



cm 

n = n GD = 1 + 0.226— p(h), 



(1) 



with p(h) being the density of air at an atmospheric 
height h. Although this index varies only between 1 and 
1.0003, this variation has important consequences, as dis- 
cussed in detail in Ref. [16]. For example, the retarded 
time t* (the time when the signal was sent out) for a given 
observer position is a multivalued function of the observer 
time t, which gives rise to "Cherenkov effect" phenomena, 
where the signal may become very short and very strong. 
The caveat in this treatment is the fact that we con- 
sider the currents to be point-like, which is only a good 
approximation far from the shower axis. The Cherenkov- 
like effects actually show up as singularities, and we ex- 
pect these singularities to disappear when we give up the 
"point-like" assumption. Nevertheless, although Ref. [l(| 
does not provide a realistic picture for all observer dis- 
tances, its results are very important as the basis of the 
much more realistic description employed in the present 
paper. 

Anyhow, the notion "point-like" has to be taken with 
care. In the point-like picture described in Ref. [l6[, we 
do not have a simple moving point-like charge, we rather 
have already transverse currents, and also the longitudi- 
nal structure is nontrivial, just all these currents are - at 
a given time - concentrated in a very small volume. But 
there must be an internal structure, and therefore it is 
natural as a next step to investigate the three-dimensional 
structure of the shower at a given time. In order to do 
so, we consider a "shower fixed" coordinate system. The 
origin O of this system is the center of the shower front. 
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Figure 1: The distribution of charged particles w(x,y,h) at a given time, as a function of the transverse coordinate x and the 
longitudinal coordinate h, for y = 0. 



We use the coordinates x and y to describe positions in 
the plane transverse to the shower axis, and h as the lon- 
gitudinal distance behind the shower front. The latter 
one is actually a hypothetical plane, which contains real 
particles only around x = y = 0, whereas for larger dis- 
tances, the fastest particles stay behind this plane. The 
situation as obtained in a realistic Monte Carlo simula- 
tion (details to be discussed later) is shown in fig. [TJ The 
distribution of charged particles shows a very sharp max- 
imum at the origin (x = y = h = 0), and falls steeply 
in transverse and longitudinal direction. We will discuss 
the functional form of this distribution later in detail, for 
the moment we only want to illustrate the fact that the 
distribution obtained from simulations shows nontrivial 
structures, concentrated in a small range in particular 
concerning the h variable. 

In the current paper, we want to take into account 
the realistic three-dimensional form (at a given time) of 
the shower, as obtained from shower simulations, still us- 
ing a realistic index of refraction. The numerical pro- 
cedures of our approach, referred to as EVA 1 (Electric 
fields, using a Variable index of refraction in Air shower 
simulations), amount to air shower simulations, analy- 
sis tools for extracting currents and shower shapes, and 
automatic fitting procedures providing smooth functions 
for all relevant shower characteristics. First results of our 
new approach have been published recently [171 . In the 
last part of the paper, we discuss important consequences 
of our approach, referred to as "geomagnetic Cherenkov 



radiation", which provides strong emissions in the GHz 
frequency domain, alone or as double peak structures in 
the frequency spectrum. 

II. TAMING SINGULARITIES 

We first repeat some elementary facts of the shower 
evolution, which have been discussed in detail in Ref. [lj| . 
We consider here showers due to a very energetic pri- 
mary particle, with an energy above 10 14 eV. Such a 
shower moves with a velocity /3c, which is very close 
to the vacuum velocity of light c. There is a constant 
creation of electrons and positrons at the shower front, 
with somewhat more electrons than positrons (electron 
excess). This is compensated by positive ions in the air, 
essentially at rest. The electrons and positrons of the 
shower scatter and lose energy, and therefore they move 
slower than the shower front, falling behind, and finally 
drop out as "slow electrons / positrons". Close to the 
shower maximum, the charge excess of the "dropping out" 
particles is compensated by the positive ions, since there 
is no current before or behind the shower. Taking all 
together we have the situation of a moving charge, mov- 
ing with the vacuum velocity of light, even though the 
electrons and positrons are moving slower, and they are 
deviated (in opposite directions) in the Earth magnetic 
field. 

Neglecting the finite dimension of the shower, referred 
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Figure 2: The dependence of the retarded time t* on the observer time t for n = 1 (dashed line), n = hod (solid line), and 
n = n groun( j ~ 1.0003 (dotted line) for inclined showers (27° and 70° and for different distances in meters of the observer from 
the impact). The reference time ts is the time of closest approach of the shower with respect to the observer. 



to as "point-like" (PL) approximation, one has a four- 
current 

j Ph {t',Z) = J(t')5 3 (x-i(t')), (2) 

with a longitudinal component due to charge excess, and 
a transverse component due to the geo-magnetic field. 
Solving Maxwell's equations, we can express the potential 
in terms of the four-current J, evaluated at the retarded 
time i*, as [H| 

4(t,2) = ai, ( 3 ) 

PL 4vr \RV\ 

with V = c d^/dt', and with R being a four- vector de- 
fined as R° = c(t - t*) and W = -Ld/d?L, where L 

is the optical path length between the source £(t*) and 
the observer. This point-like approximation is certainly 
only valid at large impact parameters (> 500 m), but even 
more importantly it will serve as a basis for more realistic 
calculations, as discussed later. It should be noted that 
in case of n > 1 and even more for n — tigd the vector 
potential shows singularities, which arise from 1/|.RV| oc 
dt* I dt and the fact that t* is a non-monotonic function 
of t, as shown in fig. [5] and discussed in detail in [16j. We 
show the realistic case n = tigd with the corresponding 
curve situated between the two limiting cases n = 1 and 
n = 1.0003. 

In general, one needs to consider the finite extension 
of the shower at a given time t' , expressed via a weight 



function w(r , r 2 ,h), where r 1 and r 2 represent the trans- 
verse distance from the shower axis, and h the longitudi- 
nal distance from the shower front, the latter one moving 
by definition with the vacuum velocity of light. Positive h 
means a position behind the shower front, and therefore 
w is non-vanishing only for positive h. The weight will fall 
off rapidly with increasing distance r — J (r 1 ) 2 + (r 2 ) 2 
from the axis. The precise form of w will be discussed in 
a later chapter. In principle one needs to convolute the 
weight w with the currents, and then compute the po- 
tential and field. Due to a translation invariance (being 
correct to a very good approximation, since the index of 
refraction varies slowly), this is the same as computing 
first the potential in PL approximation, and then per- 
forming a convolution as 

A (t,x) = Jd 2 rJdhw(f,h)A^, L (t,x^ ~ h,^ +r). (4) 

where x" and ar 1 = (a;^ 1 ,^^ 2 ) are coordinates parallel 
and transverse to the shower axis. Defining y 1 - = x + f. 
we get 

Af ) (t,x) = Jd 2 y^ Jdhwiy 1 - - x ± ,h)A p PL (t,x^ - h,y^). 

(5) 

The electric field is then obtained from the derivatives of 
A. 

One cannot simply exchange derivation and integra- 
tion, due to the presence of singularities as discussed be- 
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fore, and therefore a naive convolution of w with EpL is 
not possible: one needs a more sophisticated treatment 
of the singularities. So let us consider the most general 
case of a multi-valued function t* as a function of the 
observation time t, as sketched in fig. |3l The function 

t* I 




In the following we keep in mind that Ap L has the fol- 
lowing arguments: the time t, the longitudinal variable 
x" — hk + A, and the transverse variable y±; w has the 
arguments hk — A and y± — x±. We do not write these 
arguments explicitly, to simplify the notation. We also 
omit the arguments t, x of the vector potential. So we 
write 



A fi = d 2 v ± dXwA 



PL * 



The components of the electric field are 



= c(- 



dA° d A II 

dxW Oct 



(11) 



(12) 



t 

Figure 3: Several branches of the function of t* versus t, where 
t* is the retarded time corresponding to an observer time t . 

is composed of several branches, br„, limited by certain 
times ifc. The derivative dt* j dt becomes infinite at these 
branch endpoints, and the point-like potential becomes 
singular. This is why we refer to the tk as "critical times". 
Close to these singularities, we have 

rif* 

t* -**(**) ~ I* -**| 1/2 . and V~ |t_ * fc| ~ 1/2 - (6) 

When evaluating eq. ([5]), we have to worry about the 
critical time for a given observer position (x" — h,y^~), 
corresponding to the arguments of Apl- In other words, 
tk is a function of these variables, i.e. 



tk = tk(x" - h,y^). 



(7) 



It is useful to define a "critical h value" hk, for given t, 
via 

t = t k {J -h k ,y x ), (8) 
which allows us to write eq. ([5]) for a single branch as 

A l3 (t,x)= Jd 2 y^J "dhw^ -x ± ,h) 

(9) 



A^J-h,^). 



Using the integration variable A = hk — h, we obtain our 
master formula for the vector potential, 



h k 



A^(t,x) = d 2 y^ d\w{y^ -x^,h k -\) 



A^ L (t,J - h h + X,y^), 



which is useful because Ap^ has a singularity in A for 
A — ► 0, which does not interfere with the derivatives 
which have to be performed in order to get the fields. 



E 



dA° _ dA^_, 

dx ±l dct ' 



(13) 



Using A ph = \RV\~ 1 and eqs. (|A1IA2|) , the time 

derivative of the vector potential may be written as 



OA 1 



d 2 y_ 



</A j-K/ApL + w^ PL | 



(14) 



with w' = dw/dh, A PL = ^K l \RV\-\ K = dJ/dt'. 
In principle there is an additional term from the time 
derivative of the upper limit of integration, but this con- 
tribution vanishes due to w(0) = (see next chapter). 
Concerning the space derivative, we first compute the 
derivative with respect to the longitudinal variable. We 
find 



d 



A" 



dxW j ./,, 



j2„, / „„/ aO 



d-y± 



w' ApjdX (15) 



since the total longitudinal space derivative of A PL van- 
ishes exactly. The transverse derivatives of the scalar 
potential can be expressed in terms of the derivatives 
w l = dw/dr' 1 of the weight function w as 







8x A 



j,A°= d 2 y± dXw l A° 



1 PL- 



(16) 



The above results for the partial derivatives of the vector 
potential A^ allow us to obtain corresponding expres- 
sions for the electric field. The longitudinal electric field 
c(9"A° — d°A^ ) is given as 



= -c / d 2 y ± 



dX {w'A° PL -w'A l PL 



wAl L } 



(17) 

(10) The transverse field c(d ±l A° — d°A ±l ) can be written as 



E-' = > ld 2 y±J^ L dx{w l A P 



PL 



,u' A^-wA^ 



(18) 
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The formulas simplify considerably far from the singular- 
ity as well as at the singularity, but we keep the exact 
expressions, in order to interpolate correctly between the 
two extremes. It should be noted that the above expres- 
sion concerns a single branch, the complete field is the 
sum over all branches. 

In the present work we have derived the electric field 
directly from the Lienard-Wiechert potentials in the 
Lorentz gauge without further approximations. The dis- 
tribution of the particles in the shower front over a finite 
volume is the reason that our final result is not plagued 
with singularities in the vicinity of Cherenkov emission. 
We thus explicitly include both the near- and the far-field 
components of the radiation. In this sense it differs con- 
siderably from the calculations presented in [l4| where 
an ad-hoc frequency cut-off is introduced in the calcula- 
tions of 300 MHz, and the near-field component of the 
electric field is neglected (the Fraunhofer condition). As 
can be seen from fig. BB1 below, the data show a consid- 
erable intensity above 300 MHz. A question that arises 
in this respect is the validity of the Fraunhofer condition 
when Cherenkov effects come into play which implicitly 
is assumed in [3| ■ Often the Fraunhofer condition is for- 
mulated as a 2 sin 2 9/R < X/2n , where a is the length 
of the emission trajectory, 9 is the opening angle from 
the emission point to the observer, A the wavelength of 
the emitted signal, and R the distance from the emis- 
sion point to the observer. If there is a single point on 
a where the Cherenkov condition is fulfilled, the electric 
field will diverge at this point whereas the field is finite 
at all other points. This implies that thus the Fraunhofer 
condition is not valid. A Lorentz-invariant formulation of 
the Fraunhofer condition is a 2 sin 2 9/RV < A/2-7T, where 
the distance R is replaced by the retarded distance RV. 
Since the retarded distance vanishes at the Cherenkov 
angle this clearly shows that at this point the Fraunhofer 
condition is no longer valid for which reason we have not 
made this assumption in our approach. 



III. MONTE CARLO SIMULATIONS AND 
FITTING PROCEDURES: EVA 1.0 



The numerical evaluation of the eqs. (|17I18|) is done 
employing the EVA 1.0 package, which 

• provides the weights w, 

• provides the currents J needed to compute the po- 
tentials Ap L , 



• does the numerical integration of eqs. (|17I18[) and 
the summation over branches. 

Both the weights w and the point-like currents J are ob- 
tained from realistic Monte Carlo simulations of air show- 
ers. The EVA package consists of several elements: 



• the air shower simulation code CX-MC-GEO, in- 
cluding analysis tools to extract four-currents and 
the shape of the shower, 

• the automatic fitting procedure FITMC which al- 
lows to obtain analytical expressions for the cur- 
rents, 

• the EVA program which solves the non-trivial prob- 
lem to compute the retarded time and "the denom- 
inator" RV, for a realistic index of refraction. 

We first discuss air showers. They are considered point- 
like for the moment, as seen by a far-away observer. The 
shower is a moving point, defining a straight line trajec- 
tory, see fig. H]. One defines an "observer level" which is a 



shower position 
at time t'=0 



shower position 
at time t'=t 




trajectory 



observer G 



shower position 
at some time t' 



obs. level 



Figure 4: Air shower as seen by an observer G. The point B 
is the point of closest approach with respect to the observer 
G. The point C is the shower position at some time t' . The 
point B corresponds to the shower position at t' = is (which 
may be taken to be zero). 

plane of given altitude z with respect to the sea level. One 
defines some arbitrary point A on the trajectory. The 
corresponding projection to the observer level is named 
O (origin) and the observer position is given in terms 
of coordinates (x, y) with respect to O. The x— axis is 
the intersection of the "shower plane" OAC and the ob- 
server level. The angle between the shower trajectory 
and the vertical axis OA is referred to as inclination and 
denoted as 9. In many applications, A and O coincide: in 
this case they represent the impact point. For horizontal 
showers the two points are different. The geomagnetic 
field is specified by an angle a with respect to the ver- 
tical, and an angle ip with respect to the shower plane 
(ip = means that B points towards the shower origin). 
One can of course see it the other way round (maybe even 
more natural): for a given orientation of the geomagnetic 
field, ip defines the orientation of the shower axis. 

In the EVA framework, one has to specify the altitude 
z, the distance a = \OA\, the inclination 9, the energy E 
of the shower, and the observer coordinates x, y. And in 
addition the angles a and ip and the magnitude B of the 
geomagnetic field. 

The actual air shower simulations are done with a sim- 
ulation program called CX-MC-GEO, being part of the 
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EVA package. It is based on CONEX which has 

been developed to do air shower calculations based on a 
hybrid technique, which combines Monte Carlo simula- 
tions and numerical solutions of cascade equations. It is 
also possible to run CONEX in a pure cascade mode, and 
this is precisely what we use. This provides full Monte 
Carlo air shower simulations, using EGS4 [2(| for the elec- 
tromagnetic cascade, and the usual hadronic interaction 
models (QGSJET, EPOS, etc) to simulate hadronic in- 
teractions. 

Two features have been added to CONEX. First of 
all a magnetic field, which amounts to replacing the 
straight line trajectories of charged particles by curved 
ones. This concerns both the electromagnetic cascade 
and the hadronic one. In addition, analysis tools have 
been developed, which allow to get a complete informa- 
tion of charged particle flow in space and time. These 
features have already been developed to compute cur- 
rents in [l6j . so in particular more details about the im- 
plementation of the magnetic field can be found there 
(though we did not use the names EVA and CX-MC- 
GEO yet). We also discuss in [l6[ some details about the 
different internal coordinate systems needed to extract 
information about currents. The results shown in [l7| 
were also based on CX-MC-GEO simulations, referred to 
as CONEX-MC-GEO at the time. 

In [III, we provide several results concerning particle 
numbers and currents for different orientations of the axis 
with respect to the geomagnetic field. All the results are 
still valid, the corresponding programs did not change 
since. 

An important ingredient of our approach is the 
parametrization of the results (currents </, distributions 
w), which have been obtained from simulations in the 
form of discrete tables. This is necessary partly to per- 
form semi-analytical calculations such that numerically 
stable functions have to be dealt with without having 
huge cancellations in the results. It is especially impor- 
tant for the stable calculation of Cherenkov effects. It 
allows for the calculation of a smooth shower evolution, 
whereas when working with histogramed distributions in 
position and time, it is not possible to reconstruct a con- 
tinuous shower evolution and the artificially introduced 
sudden changes in the particle trajectories may give rise 
to spurious radio signals. 

The parametrization of Monte Carlo distributions is 
done in FITMC. This program takes the distributions 
(for currents) as obtained from the simulations in the 
form of histograms, to obtain analytical expressions, us- 
ing standard minimization procedures. FITMC creates 
actually computer code to represent the analytical func- 
tions, and this code is then executed at a later stage. The 
"basic distribution" is the so-called "electron number TV" 
(which counts the number of electrons and positrons) as 



a function of the shower time t' , which is fitted as 



AT(i') = A exp (B 



C(t' + D) 
+ E(t' + Ff 



(19) 



G{t' + H) 



As an illustration, we show here the case of a shower 
with an initial energy of 5 • 10 17 eV, an inclination 9 = 27°, 
and an azimuth angle -0 = 0°, defined with respect to the 
magnetic north pole. The angle ip refers to the origin of 
the shower. So ip = 0° thus implies that the shower moves 
from north to south. We consider the magnetic field at 
the CODALEMA site, i.e. \B\ = 47.3/jT and a = 153°, 
so the shower makes an angle of roughly 54° with the 
magnetic field. 

In figEl we plot the electron number N, as a function of 



x 10- 



:4000 -0.50E+18eV 



27.0° 0.0° 




-45 -40 -35 -30 -25 -20 -15 -10 -5 
shower time t' ((Is) 

Figure 5: The number N of electrons and positrons, as a 
function of the shower time t' for a shower with an inclination 
9 — 27° and an azimuth angle ip = 0° with respect to the mag- 
netic north pole. The full red line represents the simulation 
result, the dashed blue line is the fit. 

the shower time t', for a simulated single event, together 
with the fit curve. A thinning procedure has been applied 
(here and in the following) to obtain the shown simulation 
results. The time t' = refers to the point of closest 
approach with respect to an observer at x = 0, y = 500 m, 
z = 140 m. We suppose a = (so the shower hits the 
ground at x = y = 0). 

The magnitudes of the components J M of the currents 
have a similar t' dependence as N(t'). Therefore we 
parametrize the ratios J^ 1 /{Nec), with N being the elec- 
tron number, e the elementary charge, and c the velocity 
of light. We use the following parametric form: 



N(t')ec 



= A + B(x + C) + D(x + E) 2 + F(x + G) 3 . (20) 



In fig. |51 we plot the longitudinal current component J z 
(divided by Nec), as a function of the shower time t' , 
for a simulated single event, together with the fit curve. 
At early times - far away from the shower maximum - 
there are of course large statistical fluctuations. But since 
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Figure 6: The longitudinal current component J z , divided by 
Nec, as a function of the shower time t'. The full red line 
represents the simulation result, the dashed blue line is the 
fit. 



N(t') is very small here, this region does not contribute 
to the pulse. In fig. [JJ we plot the transverse current 
components J x and J v , (divided by Nec), as a function of 
the shower time t' , for a simulated single event, together 
with the fit curves. 

The EVA program uses these analytical fit functions 
for the current components, 



e ■ c, 



(21) 



to compute the vector potential. The currents have to 
be evaluated at t' = t* , the retarded time. The cen- 
tral part of EVA is actually the determination of the re- 
tarded time t*(t, x) for a given observer position. This is 
quite involved - in case of a realistic index of refraction 
- and described in detail in gain without referring 

to EVA, but these are exactly the same programs being 
used). A results of such a calculation is shown in fig. [5TJ 



27.0 0.0 




-45 -40 -35 -30 -25 -20 -15 -10 -5 
shower time t' (p.s) 



Figure 7: The transverse current components J x (lower lines) 
and J y (upper lines), divided by Nec, as a function of the 
shower time t' . The full red lines represents the simulation 
result, the dashed blue lines are the fits. 



A new feature compared to [16| - and most relevant 
for this paper - is the possibility to obtain information 
about the shape of the shower via the weight function 
w. The weight function w is not perfectly cylindrically 
symmetric, due to the geo-magnetic field but also due to 
statistical fluctuations, since we are considering individ- 
ual Monte Carlo events. However, in this paper we will 
neglect these tiny deviation from symmetry, and consider 
a weight function w(r, h) depending only on the two vari- 
ables r and h, related to the general weight function as 



w(r, h) = 2nr w{r, h) . 



(22) 



The lateral coordinate r measures the transverse distance 
from the shower axis, the longitudinal coordinate h is 
meant to be the distance behind the shower front. This 
front is a hypothetical plane moving parallel to the shower 
axis with the velocity of light c, such that all the particles 
are behind this front, expressed by a positive value of h. 
We will express the weight function as 



w(r, h) = w\ (r) W2 (r, h) , 



(23) 



with J drw\{r) = 1, and with J dhw%{r, h) = 1 for all 
values of r. 

We use again CX-MC-GEO to obtain w, then FITMC 
to obtain an analytical function, which is later used in the 
EVA program to compute the fields, based on the formu- 
las described in the preceding chapter. All the simulation 
results shown in the following are based on the the same 
shower, mentioned earlier when discussing currents. 

We first investigate the radial distribution wi(r). In 
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r(m) 



Figure 8: The radial distribution wi(r). The thin yellow lines 
correspond to different times, the points represent an average, 
and the thick red line corresponds to a fit (see text). 

fig. [SI we show the radial distribution as obtained from 
the Monte Carlo simulation. The thin lines correspond to 
different times t' , between —25 /is and —5 /is. The points 
represent an average over all times, and also averaged 
over r-bins. Since the time dependence is quite small, we 
will use the radial distribution at the shower maximum 
'•max as time-independent distribution w±(r). The thick 
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red line corresponds to a fit to the Monte Carlo data, 
using the form 



wx(r) 



r(4.5 - s) 
r(s)r(4.5-2s) \r 



-4.5 



(24) 

with fit parameters ro and s (providing an excellent fit). 

Knowing wx(r), we now investigate how far the par- 
ticles are moving behind the shower front, expressed in 
terms of the longitudinal distance h, for a given trans- 
verse distance r. From the above simulation, we obtain 
easily the mean distance h at a given r. We find a per- 
fectly linear time dependence, of the form 



h = h 



front 



cA/3t', 



(25) 



where A/3 can be obtained from fitting time dependence 
at different distances r, the result is shown in fig. |H] 
as solid line. The quantity A/3 represents the veloc- 




r(m) 



Figure 9: The longitudinal velocity difference A/3 versus r. 
We show the results for realistic simulations (thick red solid 
line) and for 7 = 60 (green dotted line). Also shown: the 
value 1 — l/ng roU nd (blue dashed line). 

ity difference (in units of c) with respect to the the 
shower front, which itself moves with velocity c. So 
the velocity of the "average position" of the shower is 
1 — A/3. Also shown in fig. [5J as dashed line, is the value 
1 — l/n groun d, corresponding to the velocity of light in 
air with n gro und = 1-0003. And we also plot as dotted 
line the A/3 obtained from 7 = 60, corresponding to the 
average electron energy. The simulated curve (thick full 
line) is considerably below this dashed and the dotted 
curves, which means that the velocity of the average po- 
sitions is larger than c/n groun d, it is also larger than the 
velocity of the average electron. The simulated velocity is 
even (slightly) larger than c. This is due to the fact that 
matter is moving on the average from inside (small r) to 
outside (large r), and the average h decreases with de- 
creasing distance r. But the effect is small, the deviation 
of the shower velocity from c is less than 1/1000. We 
will ignore the small time dependence for the moment, 




r(m) 



Figure 10: The mean value h for given values of the lateral 
distance r. 
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Figure 11: The parameters Hi, H2, H3, H4, and if as a 
function of the lateral distance: Hi (full line) , H2 (dashed 
line) , H3 (dotted line), H4 (dashed-dotted line) , K (wide- 
dotted line) , 



and consider in the following quantities at t max . To get 
some idea about the typical scales of the /i-distribution 
W2(r,h), for a given value of r, we determine the mean 
value h, as shown in fig. 1101 The mean value h is almost 
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Figure 12: The parameters H$ and K' as a function of the 
lateral distance: K' (full line) , H5 (dashed line) 
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a linear function of the distance r, and for r = 100 m we 
get an average h of roughly 10 m. 

The W2 distribution is obtained by fitting Monte Carlo 
data in a range h between zero and 5 h, for given r. We 
use 



w 2 (r, h) 



MGD 
IGD 



III 



(r, h) for r > ro 
(r, h) for r < tq ' 



(26) 



with uj^ 100 being a "modified gamma distribution" of the 
form 



MGD 



(r,h) 



H{r, h) G(r, h) 
N{rj ' 



(27) 



with 



H(r,h) = QiH.-h) [ 2 ( A) - (A) j +6^-^), 

(28) 



and 

G(r,fc) 



e(H 3 ~h) (h K - 1 e- h ' H ^ , (29) 

+ 0(/l - ff 3 ) (H; K - 1 e -H 3 /H 2e -(h-H 3 )/H^ > 



with AT being a normalization constant such that 
J dh W2(r, h) 

gamma distribution" of the form 



1. The function w 2 GI) is an "inverse 



„IGD 



(r,h) 



K 



,-K'-l. 



-H B /h 



We use ro = 20m. The r-dependence is hidden in the pa- 
rameters Hi, H2, H3, H4, H$, K, and K' . In figs. [TTIand 
[T2l we plot the parameters, as obtained from fitting the 
Monte Carlo data. All parameters grow with increasing 
distance r. Whereas H2 seems to saturate, all the other 
parameters grow roughly linearly with r. With these pa- 
rameters, we get good fits for h values up to five times the 
mean. In figs. [T3]and[131 we show the fits of w 2 together 
with Monte Carlo simulation results for different times. 
In fig. [TSJ we show the fitted W2 curves for three differ- 
ent values of r, conveniently plotted as hwz versus h/h, 
where one clearly sees the evolution of the shape with r. 

The reason to switch between w^ 100 and u;!, 00 & t 
ro = 20 m becomes clear from figures [T3] through [15j 
From figure [15] it can be seen clearly that the particle 
distribution as obtained from the Monte-Carlo simula- 
tions behaves quite differently close to the shower axis as 
compared to the distribution at large distances. This dif- 
ferent behavior requires the use of different fit-functions 
in both regimes. At large distances from the core, the 
parametrization of w^ GD reproduces the MC result ac- 
curately. At small distances, it is important to have a 
smooth parametrization without jumps in the first deriva- 
tives, which is the case when using w 2 GI} - 
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Figure 13: The distribution W2(r,h) for r = 5m. The full 
black line represents the fit, the dotted lines are simulation 
results for different times. 
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Figure 14: The distribution W2 (r, h) for r = 128 m. The full 
black line represents the fit, the dotted lines are simulation 
results for different times. 



The above fit function w 2 D leads to a delta-peak at 
r = 0. To still obtain numerical stability, a cut-off for 




Figure 15: The distribution W2(r,h) for r = 5m (full line), 
r = 30m (dashed line), and r = 180m (dotted line). 
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the values K' and H$ is introduced such that the width 
of u>2 GD is 1 mm. Since most of the particles are located 
within r = Sx 1 - = 1 m from the shower axis, the path dif- 
ference between signals emitted at this distance on both 
sides of the shower axis acts as the important length scale 
in this regime. We estimate this path difference SR. for 
a constant index of refraction equal to n = 1.0003 : we 
have SR ~ -^rrSx 1 - ~ \J n 2 /3 2 — lSx 1 - « 3 cm. Here we 
use that at the Cherenkov time (critical time for h = Q), 

we have R° — n/3x", and x'j = ^n 2 /3 2 — Ix^r [2l|. So a 
cut-off of wl GB « 1 mm should give stable results. This 
has been tested numerically. 



IV. TIME SIGNALS 
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Figure 17: The x component of the charge excess contribution 
to the electric field as a function of the observer time t in 
ns, for an observer distance of 112 m (upper panel) and 448 
m (lower panel) We compare different options for the index 
of refraction n, namely n — l(left), n — uqd (middle), and 
n = 1.0003(right). 



and n — 1.0003. One can clearly see big differences be- 
tween the three scenarios, up to a factor of ten in width 
and magnitude. We also see, even in the realistic case 
(n = «gd), the appearance of "Cherenkov-like" behavior, 
with very sharp peaks. In figs. [TS] and 1191 we consider 
a more inclined shower (70°), for two different observer 
positions: 292 and 1170 meters to the south of the im- 
pact point. The differences between the realistic case 



Figure 16: The y component of the geomagnetic contribution 
to the electric field as a function of the observer time t in 
ns, for an observer distance of 112 m (upper panel) and 448 
m (lower panel) We compare different options for the index 
of refraction n, namely n = l(left), n = nGD(middle), and 
n = 1.0003(right). 



As already said, the eqs. (|17ll8p are evaluated em- 
ploying the EVA 1.0 package, which provides the weights 
w, the currents J, the denominators RV, and the inte- 
gration procedures, as discussed in the previous chapter. 
We first consider the same "reference shower" (initial en- 
ergy of 5T0 17 eV, inclination 27°) already discussed there. 
We will distinguish between the geomagnetic contribution 
(caused by the currents due to the geomagnetic field) and 
the contributions due to charge excess. In figs. [TH] and 
[T71 we show the results for the two contributions, for two 
different observer positions: 112 and 448 meters to the 
south of the impact point. We compare the realistic sce- 
nario (n = «gd) with the two "limiting cases" n = 1 
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Figure 19: Same as fig 1171 but here we consider a more in- 
clined shower. 

(n = ugb) and the two "limiting cases" is even bigger: 
more than a factor of 100 in width and magnitude! 

V. GEOMAGNETIC CHERENKOV RADIATION 

As shown in the last chapter, a realistic treatment of 
the index of refraction in the atmosphere seems to be very 
crucial for the forms of the electromagnetic pulses. Can 
this be seen in experiments? What exactly should one 
look for? 

To answer these questions we will discuss in the follow- 
ing frequency spectra. As shown in chapter II, the fields 
are sums of terms of the form (up to factors) 

J dV pancake x currents x (^~^J > (30) 

where "currents" and "pancake" refer to respectively the 
pointlike currents and the current distributions in the 
pancake, or its derivatives. The quantity dV is a pan- 
cake volume element. The currents and the "Cherenkov 
term" dt*/dt are taken at the retarded time t' — t* , for 
a given observer time and position. Let us consider the 
evolution of an air shower in time t' . The currents are 
essentially proportional to the electron number N e (t') of 
the shower, the so-called "profile". We define t* p to be 
an emission time (retarded time) corresponding to the 
profile maximum, also referred to as shower maximum. 
Another important quantity is the Cherenkov time t£», 
corresponding to the time where dt*/dt becomes singu- 
lar. 

The electric field contains actually terms governed by 
the derivatives of the currents, and therefore by the 



derivative of the profile. We consider therefore the ex- 
pression "shower maximum" to represent the actual max- 
imum of the profile or of its derivative. 

A strong signal is expected when the two times t* p and 
t* c coincide. Such a situation is shown in fig. 1201 where 




frequency (MHz) 



Figure 20: Fourier transform of geomagnetic component of 
the 70 degrees inclined shower observed at 1170 meters from 
the shower core. We plot the modulus of the Fourier trans- 
form. 

we plot the Fourier spectrum for the geomagnetic com- 
ponent of the electromagnetic field for the 70 degrees in- 
clined shower discussed in the previous chapter, with an 
initial energy of 5 • 10 17 eV, and an observer positioned at 
a distance of d = 1170 m to the east of the impact point 
of the shower, corresponding to an impact parameter b 
of around 400 m. At this distance the shower maximum 
occurs at the Cherenkov time for a realistic index of re- 
fraction. The realistic case (n = ugd) contains very high 
frequency components up to several GHz as one would ex- 
pect from the sharp peak in figure 18. The two limiting 
cases peak at lower frequencies below 100 MHz. 

In the following, we discuss some very interesting fea- 
tures by taking the example of a 60 degrees inclined 
shower with an initial energy of 10 17 eV, moving from 
west to east, in a magnetic field of strength 24.3/xT and 
an inclination a of 54° (Auger site). The observer is po- 
sitioned to the east of the impact point. We will use the 
impact parameter rather than the horizontal distance (as 
in the examples before) to characterize the observer po- 
sition. 

In fig. [2TJ we plot the shower profile N e as a function of 
the retarded time £*, together the retarded times t* as a 
function of the observer time t, for three different choices 
of the impact parameter. For large values of b (above 
285m), like the case of 300 meters (magenta curve), there 
is no Cherenkov time, the function t*(t) is single val- 
ued, the derivative is always finite. We have "normal" 
emission, coming from around the the maximum of the 
profile corresponding to t' = tp, see fig. [33] The form 
of the time signal is determined by the profile, we expect 
maximum frequencies around few hundred MHz, as con- 
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Figure 21: The shower profile as a function of t* (black line) 
and the retarded times t* as a function of the observer time 
t, relative to the time of closest approach ts (red, blue, and 
magenta curves). The "Cherenkov points" correspond to the 
Cherenkov times (where dt* jdt is singular). 



firmed by the calculation shown in fig. 1251 At impact 
parameters smaller than 285 meters, the function t*(t) 
starts to become double valued, so we we start observing 
a Cherenkov time. At 250 meters, the Cherenkov time co- 
incides with the shower maximum, we have Cherenkov 
emissions from around the shower maximum. This means 
that due to dt* jdt = oo, the emissions from a broad re- 
gion around the maximum will be "compressed" and ar- 
rive almost at the same time at the observer, as sketched 
in fig. [23l This leads to a strong and very short signal. 
Since the singularity is integrated over, as explained in 
chapter II, the actual width of the same signal is deter- 
mined by the current distributions in the pancake, and 
we expect frequencies around a GHz, as confirmed by the 
calculation shown in fig. [25] 

If the observer is even closer to the shower, for exam- 
ple at an impact parameter of 180 meters, we still have 
a Cherenkov time, but this time is now significantly later 
than the shower maximum time (see the dot on the red 




O 




"o 



Figure 23: The observer O receives "Cherenkov emission" 
from around the shower maximum. 




Figure 24: The observer O receives both "normal emission" 
from around the shower maximum and "Cherenkov emission" 
from later times. 



curve in fig. |2"T|) . Here we may have a very interest- 
ing situation: the observer may receive "normal" emis- 
sion from around the shower maximum, but at the same 
time he may receive a significant contribution from much 
later, around the Cherenkov time, which again due to the 
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Figure 22: The observer O receives "normal emission" from 
around the shower maximum. 



Figure 25: Flux densities for radio emission from a 10 17 eV 
energy shower at 60° zenith angle for impact parameters of 
180, 250, and 300 meters. 
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Cherenkov effect (signal compression) will be relatively 
strong and short (high frequency, order of GHz). This 
situation is sketched in fig. [53] The calculations in fig. 
1251 show (as expected) two distinct peaks, one at small 
frequencies due to the normal emission from the shower 
maximum, and a second peak at high frequencies due to 
Cherenkov emission at much later times. 

So our approach predicts not only high frequency com- 
ponents due to the geomagnetic Cherenkov effect, but in 
addition a double peak structure which reflects the simul- 
taneous reception of signals from very different positions 
of the shower: "normal" emissions from around the max- 
imum, and Cherenkov emission from much later times. 



VI. COMPARING TO DATA 

This geomagnetic Cherenkov radiation might have 
been observed by the ANITA-collaboration |22||, where 
pulses have been measured in the 200-1200 MHz band. 
Furthermore, since these high frequency components oc- 
cur only at the Cherenkov distance, upon applying a high- 
pass filter a clear Cherenkov ring should become visible 
in the LDF. The radius of this ring contains direct in- 
formation of the shower maximum and thus the chemical 
composition of the original cosmic ray. New experiments 
at the Pierre Auger observatory [1, [f|, and LOFAR [HJ 
should be able to measure the LDF in more detail, where 
first hints of "Cherenkov- like" effects might have been 
observed (1[2^|. 

The key result of our present work is the prediction 
of a sizable power emitted at higher frequencies, and a 
possible double peak structure with one peak at high fre- 
quencies. There exist only few observations where the 
spectrum over a large frequency range has been mea- 
sured. A good example was published recently by the 
ANITA collaboration [25||, showing the summed power of 
two cosmic-ray events for the range of 300-900 MHz. 

In this measurement no indication is given of the arrival 
direction and the energy of the initiating cosmic ray, only 
that it most probably came from a relatively large zenith 
angle. The azimuth angle is unknown. Therefore also the 
air density along the path of the air shower is unknown, as 
well as its orientation with respect to the magnetic field. 
All this makes any quantitative comparison impossible. 
To get at least a qualitative understanding, we compare 
the data with the result of a simulation for a cosmic-ray 
at a zenith angle of 60° , moving from west to east, in 
a magnetic field corresponding to the Auger site, with 
an observer east to the impact point, for various impact 
parameters b - the same situation as discussed in the 
previous chapter (changing the energy and the arrival 
direction of the cosmic ray will not change the qualitative 
discussion). 

In fig. 1261 we compare the data with our simulation 




200 300 400 500 600 700 800 900 1000 

KMHz) 

Figure 26: The predicted flux densities for radio emission 
from a 10 17 eV energy shower at 60° zenith angle for various 
impact parameters b are compared to the data for the sum of 
two events as measured by the ANITA balloon mission (25| . 
where the data are taken from fig. 3 of that publication. 

results. We show blue curves corresponding to 350 - 225 
m, from bottom to top for the leftmost value. The red 
curves refer to 200 - 180m, from top to bottom. 

From the discussion of the last chapter, we easily un- 
derstand the different theoretical curves: for large im- 
pact parameters 350, 325, 300m) we have the situation 
corresponding to fig. 1221 normal emission from around 
the shower maximum dominates. For impact param- 
eters around 250 meters, we have Cherenkov emission 
from around the shower maximum, as in fig. 1231 we 
get strong signals at large frequencies (GHz). Then fi- 
nally below 200 meters, we have the situation sketched 
in fig. [2U a double peak structure due to simultane- 
ously arriving signals from very different positions of the 
shower: "normal" emissions from around the maximum, 
and Cherenkov emission from later times. 

Although energy and inclination of the measured show- 
ers are unknown, it is nevertheless clear that the data 
show a structure similar to the transition region towards 
a double peak behavior, as predicted in our calculations 
shown in fig. [2H 

VII. SUMMARY 

We presented a realistic calculation of coherent electro- 
magnetic radiation from air showers initiated by ultra- 
high energy cosmic rays. The underlying currents are ob- 
tained from three-dimensional Monte Carlo simulations 
of air showers in a realistic geo-magnetic field. We take 
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into account the correct shape of the particle distribution 
in a shower at a given time. The numerical procedures - 
simulations, fitting procedures, convolutions, referred to 
as EVA 1.0 - have been discussed. We showed the im- 
portance of a correct treatment of the index of refraction 
in air, given by the law of Gladstone and Dale: using 
the correct index of refraction tiqd gives very different 
results compared to a simplified treatment using a con- 
stant index, with differences in width and magnitude up 
to a factor of 100. The new treatment leads in particu- 
lar to important emission at high frequencies (GHz). In 
certain cases, double peak structures are predicted, due 
to signals arriving simultaneously from different positions 
of the shower: "normal" emissions from around the max- 
imum, and Cherenkov emission from later times. 



Appendix A: Some derivatives 

Theorem: The quantity h k is a function of t and 
x " (the transverse coordinates are not considered here). 
Its derivatives are 



d , 

-—h k = -1, -r-uhk = 1 



dct 



dxW 



(Al) 



(A2) 



The time derivatives of ct* and RV vanish: 

= *r V = o. 
act act 

Proof: 

In the following, we do not consider explicitly the trans- 
verse coordinates (to be considered constant). The vari- 
ables of interest are the time ct = x° and the longitudinal 
coordinate x" = z. We use for any function f(t, z) the 
notation d a f = df/dct, d z f = -df/dz. For the "total" 
derivatives, we use d°f = df/dct, d z f = —df/dz. 

We consider for given fixed h the retarded time t*(t, z — 
h) . The retarded time corresponding to a critical time is 
given as 

t* k = t*(t k ,z-h), (A3) 
with tk = tk(z — h). So we have t* k =t* k (z — h). We have 

d z t* k = d°ct*d z t k +d z t*. (A4) 

In the following we make extensively use of definitions 
and relations from reference [Hj]. For t — t k , we have 

RV = c(t k -t* k )+R j V j \ t . =tt =0. (A5) 
We compute the derivative d z : 

= d z ct k - d z ct* k + {g z 3 - Vjd'ctDV 1 (A6) 
= d'ct k - (1 - V j V j )d z ct* k + V z (A7) 
= d z ct k - VVd z ct* k + V z (A8) 
= d z ct k - VV(d°ct* d z ct k + d z ct*) + V z , (A9) 



which leads to (using V z — V") 

-V II +VVd z ct* 



d z ct k 



l-VVd°ct* 



Using 



we get 



d a ct* 



R a 

W 



d z ct k = - 



yii - VVRW/EV 
1-VVR°/RV ' 



(A10) 



(All) 



(A12) 



The time derivative of h k as obtained from its definition 
is 



d°h k = {d'ctk)- 1 



which gives 



d°h k = - 



RV - VV R° 



RVVW - VVRW 
Using RV = 0, and M w L = R°, we find 

dPh k = -1. 
The other derivative of h k is trivial: 
d z h k = -1. 



(A13) 



(A14) 



(A15) 



(A16) 



The potential A, its denominator RV, and also the ar- 
gument ct* of its numerator J are evaluated at the parallel 
coordinate x" — h k + A, so the total time derivatives are 



explicitly 



d° =d° + d°h k d z , 



d° = d°-d z . 



(A17) 



(A18) 



We get 



d°ct* = (d° - d z )ct* 



R° -R z 
RV 



and with 

d a RV = V a - VVd a ct* 
(eq. (22) from [l6j]), we obtain 

d°RV = 0. 



0, (A19) 



(A20) 



(A21) 
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